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Abstract 

We present all the possible solutions of a Josephson junction with bias current 
and magnetic field with both inline and overlap geometry, and examine their 
stability. We follow the bifurcation of new solutions as we increase the junction 
length. The analytical results in terms of elliptic functions in the case of inline 
geometry, are in agreement with the numerical calculations and explain the 
strong hysteretic phenomena typically seen in the calculation of the maximum 
tunneling current. This suggests a different experimental approach based on 
the use, instead of the external magnetic field the modulus of the elliptic 
function or the related quantity the total magnetic fiux to avoid hysteretic 
behavior and unfold the overlapping Imax[H) curves. 



I. INTRODUCTION 



The static properties of a narrow Josephson junction are well characterized by the static 
sine-Gordon equation [|l|. They are experimentally measured by the maximum tunneling 
current Imax as a function of the external field H. This is an important and useful mea- 
surement since it is required not only for the characterization of the junction properties 
but also to tailor a device with the desired maximum current. In contrast to its simple 
form, the static sine-Gordon differential equation problem poses several mathematical and 
computational challenges. The complete analysis of all of its solutions is hard due to vari- 
ous interesting properties (nonlinearity, non definiteness, periodicity, boundary conditions of 
Newmann type) inherent to the sine-Gordon problem and the determination of Imax either 
theoretically, numerically or experimentally is difficult. 

Several studies analyzing the dynamic and static stability of fiuxons in the sine-Gordon 
equation have appeared in the literature in the past two decades with being the most 
representative of them. All these studies combine theoretical and numerical analysis and 
they mainly address the case where there is no external current or magnetic field applied 
on the junction. None of the above studies is comprehensive enough as far as exploiting all 
solutions and studying the affect of all physical and geometrical parameters of the problem. 

The main objective of this study is fourfold. 

• To analytically express all static solutions of 1-dimensional narrow Josephson junctions 
in a way that will allow us to examine their stability properties and their evolution 
with respect to the size of the junction, and the applied magnetic field and current. 

• To explain the hysteretic behavior and if possible to find the important physical pa- 
rameters that unravel the hysterisis. 

• To build a numerical simulation framework that will allow us to verify some of our 
theoretical results and show that they apply to more complicated Josephson junction 
configurations. 

• To propose an experimental procedure that will enable us to examine the properties 
of superconducting devices in a more accurate way. This approach is particularly 
useful for the analysis of devices that deviate from the standard mathematical model 
currently used i.e. junctions with impurities, inhomogeneities, ... 

We should mark here that, at the present, a complete theoretical analysis of such devices 
is not feasible while numerical simulations, based on state-of-the-art software packages , 
usually fail to capture all the important features. This is mainly due to the difficulty to 
track the continuation branches and to deal with the bifurcation points involved. It is worth 
to mark here that the effects of the above mentioned difficulties are clearly seen, even in a 
relatively simple case, by the fact that it was only very recently [^] realized that the critical 
value of the bias current corresponds not to a termination point, as conjectured for many 
years, but to a turning point in the bifurcation diagram. 

In addition to the above, the experimental analysis and the computer simulation and 
analysis of superconducting devices modeled by the sine-Gordon equation require an initial 
guess that is reasonably close to the desired solution. The selection of such initial guesses 
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significantly affects the effectiveness of the various continuation methods needed to determine 
Imax- These guesses will be extensions of the various 1-D solutions obtained here. Therefore 
the results of the present study are expected to be fully utilized in effectively analyzing two 
dimensional window Josephson junctions. 

The behavior of a Josephson junction is determined by the Josephson characteristic 
length Aj which depends both on material and geometry properties. For a short junction of 
length / -C A J the Imax vs H pattern (shown in Fig. |I|a for normalized length to = ^ = 1.0) 
presents the usual Fraunhoffer pattern 

T _ £0_ 

max <j) 5 

where $ = Hid is the applied flux, $o = ^ the flux quantum and d the magnetic thickness 
III]. Each of the lobes in the diagram can be labeled by the pair of integers (n, n + 1) where 
at one end the magnetic field corresponds to exactly n fiuxons (i.e. fiux= n$o) and at the 
other end n + 1 fiuxons. For the case of a long junction the problem was solved by Owen 
and Scalapino 0, and there, the different lobes overlap (as shown in Fig. [l|b for w = 10). 

It should be remarked that the sine-Gordon equation, due to its nonlinearity and pe- 
riodicity with several equilibrium points, has a multiplicity of solutions as shown by the 
overlapping lobes (in Fig. ||b) and the existence of several unstable branches which play an 
important role in the hysteretic behavior as you vary the external magnetic field. The unsta- 
ble branches are of interest too because one can stabilize them, by introducing small defects, 
and therefore they should lead to observable maximum current. As we will see later in the 
discussion a defect can modify significantly the relative amplitude of the different lobes in 
the Imax vs H curve. The unstable branches can be partially traced experimentally if we 
perform a quasistatic scanning of the magnetic field. 

The rest of the paper is organized as following. In section 2 we present the mathematical 
problem, give the explicit analytical solutions using elliptic functions and sketch the stability 
analysis. In section 3 we study the solutions in the particular case of zero magnetic field 
11 = 0, and in section 4 the analytic solutions for zero current / = 0. We study both 
theoretically and numerically their stability in section 5 and we calculate the maximum 
tunneling current. In section 6 we briefly propose an experimental procedure which utilizes 
our numerical procedure. We summarize our results in the last section. 



II. INLINE GEOMETRY AND STABILITY OF SOLUTIONS 

The case of inline geometry is a one dimensional problem, even for a wide junction, and 
one can obtain analytic solutions §]. Furthermore one can easily check their stability by 
using linear perturbations. The particular cases of zero external magnetic field and inline 
bias current also reduce the overlap boundary conditions to inline ones. In the 1-dimensional 
case we have to solve the following problem 

with the inline boundary condition 
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J rr 

x=±w/2 = ±- + ii = 7±, (2) 



dx 
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where w is the normahzed junction length and / is the current hne density. 

Eq. has a static solution ^q{x), implicitly expressed using eUiptic functions as 

sin q{x) = —2^/m sn {x + XQ\m) dn{x + XQ\m), (3) 
cos $0(2;) = 2m sn^ (x + Xolm) — 1, (4) 

2\/rn cn {x + XQ\m) , (5) 



dx 

where the modulus m determines the period of the cn elliptic function (equal to 4:K[m)) 
and the arbitrary constant xq the phase at the center of the junction. They are determined 
by the boundary conditions (^. Introducing (^) into (Q) we get 

cn{xQ + —\m) = (6) 

2y/m cn{xo — -^|"^) = 7_.. (7) 

A useful quantity to classify the solutions is the fluxon content or the magnetic flux in 
units of the quantum of flux, defined by 

At specific values of H, Nj takes integer values, so that the flux is that of an integer number 
of fiuxons. 

To check the stabihty of the static solution we consider small perturbations on 

the static solution $o(^) in the form 

<l>{x,t) = %{x) + U{x,t) (8) 

and linearize Eq. (|I|) with respect to U{x,t), to obtain 

- Utt + U^^ = cos ^o{x)-U. (9) 

There is no loss of generality if we consider specific perturbations in the form 

U{x,t) = X{x)e''. (10) 

This way we obtain from Eq. (^ the eigenvalue equation 

- X" + cos ^o{x) ■ X = X- X, (11) 

under the boundary conditions 
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X'\x=±w/2 — 0. 



(12) 



In Eq. ( pTl) A = —s^. It is seen from Eq. (0) that if the eigenvalue problem ([Tl|)-(p!2D 
has a negative eigenvalue the static solution ^o{x) is unstable. If all the eigenvalues are 
positive ^q{x) is stable, while the case A = corresponds to neutral stability and defines the 
boundary of stability. In the following we consider the two special cases 7+ = 7_ (/ = 0) 
and 7+ = — 7_ {H = 0) separately, since the associated boundary conditions are easy to 
handle and the stability analysis is simplified. 

To get a feeling concerning the possible solutions, we plot from the boundary conditions 
(0) and (^0) the constant H and I contours in the plane of the m and Xq parameters in 
Fig. ^ We give different plots for m < 1 and m > 1. The lines labeled with "0" correspond 
to if = (or 1 = for I contours) and in both cases their network encloses areas with 
a single maximum (denoted by "+") or minimum ("— ") inside. Notice that there are two 
types of curves on which H = (or / = 0) as summarized in Table |. 

The curved lines in Fig. ^ correspond to the solutions of the first line in the table which we 
call fixed Xq solutions, in the sense that the shift is a fixed fraction of the period. From (^^ 
we see that the physical quantities / and H are periodic functions of Xq with period AK^m). 
There are also solutions that have a fixed m = m* and arbitrary xq, and correspond to the 
vertical lines in Fig. § (remark that contour fitting can be distorted when two equicurrent 
lines cross). Similar results hold for the constant / contours. For m > 1 at the vertical 1 = 
curves we have an integer number of fiuxons in the junction. Thus on the lines through c (a) 
we have Nf = 1 (2) correspondingly. In the current contours the points / and b correspond 
to peaks in the Imax (see section 4.4). In the following we will focus our attention on the 
solutions where either H 01 I . 



III. NO MAGNETIC FIELD CASE (H=0) 



In the absence of external magnetic field we have 7+ = — 7_ = | and Eqs. (|,^ reduce 



2 

to 



2\/mcn (xo + ^|m^ = ^, (13) 





w , \ 




m 


2' J 



2- 



which determine the parameters m and Xq that characterize the periodicity and phase shift 
of the static solutions. There are two different classes of solutions due to the antisymmetry 
of the boundary conditions, which can be satisfied for different / either by fixing xq or m. 



A. Fixed xq solutions 

It is seen from Eqs. (0,0) and the symmetry of the elliptic functions that for positive 
/ one choice for xq in (|T^,|T^) is xq = —K{m) {K{m) is the elliptic integral of the first kind), 
so that Xq is fixed to | of the period of the elliptic function. It is in that sense that we 
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call them fixed xq solutions. Strictly xq is not a constant independent of m since K{m) is a 
function of m. Thus we have (see |^) 



I I — : sn[^\m] , ^ 

i^V^d-'") 0S"'<1. (15) 

C?l^ fx I TTT-l 

cos $o(a;) = 2 m , ^, ; - 1, < m < 1, (16) 
an^[x\m) 

i CTl { X TTl ) 

sin$o(a;) = 2 Jm (1 - m) , ', < m < 1. (17) 

Another possibility is xq = K{m) {xq is shifted by half of a period) for the case that 
2K{m) < ^ < 4:K{m), or more generally when sn{^\m) < 0, since we are limiting ourselves 
to / > 0. This means that every time w increases by 27i we introduce two extra solutions. To 
put it in other words the function I{m), in (|T5|) , is highly oscillating for large w. We do not 
need to consider the case that m > 1 since in that case we cannot satisfy the antisymmetric 
boundary conditions for the external current. 

In Fig. ^ we present three plots of J vs m for ti? = We see that for small 

w < 271 there is, as expected, only one lobe for xq = —K{m), and as we will see later only 
the part to the right of the maximum will correspond to stable solutions, while the peak 
corresponds to the maximum current for zero magnetic field. For w = ^ we have an extra 
lobe, with the left one at xq = K{m) and the right at xq = —K{m). For w = 10 (dashed 
line in Fig. ^d) the right lobe has a maximum within 10"'^ of m = 1 and corresponds to 
Xq = K{m), while the left to Xq = —K{rn). The jump in path along two different curves is 
necessary because of the restriction of positive current J. Of course the curve is symmetric 
about the J = line. The right lobe for w = 10 is shown in the inset of Fig. in expanded 
form (solid line) with a different scale for m. The part that is of experimental interest is the 
last lobe near m = 1 to the right of the maximum. The two extremal m values correspond 
to trivial solutions ^o{x) = vr (m = 0) and ^o{x) = (m = 1), the first of which is clearly 
unstable (pendulum analogy) and the second is stable. For currents above zero at w = 10 
there are four possible solutions for a given J < /*, where /* is the maximum of the lowest 
lobe. 

Because the last lobe for large w is very steep it is useful to give some analytic formulas 
valid near m = 1 and for the maximum point. By using asymptotic formulas and assuming 
that nil = 1 — m ~ e^, where e is a small parameter, we obtain the value of mi where / is 
a maximum as 

4 



^max 



^ sinh^ f ■ 
The result is consistent with our scaling assumption with 

1 



sinhf 



sm e. 



Thus care is required when simplifying the analytic formulas in (see p. 574). The 
corresponding maximum current is 
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^rff' ^^^^ 



so that for large junction length w it approaches exponentially the infinite length limit. To 
the right of the maximum the relation between /(mi) is 



3/2tanhf 



^ sechf ■ 

From the previous discussion we see that as we increase the junction length w we obtain 
more solutions. In Fig. ^ we give as a function of w the range of m values for each type of 
solution and the separating lines. These values were determined by solving numerically. 
We remark that consecutive pairs of regions of solutions correspond to different xq (i.e. 
different lobes of Fig. |^). We see that when w increases by 2tt a new pair of solutions is 
introduced. Thus near w = 10 we have four solutions labeled u, al, ar and a^, the first two 
corresponding to Xq = K{m) and the last two to Xq = —K{m). For w ^ oo there is an 
infinite number of solutions and many of the dividing lines coalesce at m = 1. The stability 



is checked by looking at the eigenvalues of the linearized problem in (pUf) . We see that 
already for w = 10 only a small range near m = 1 gives stable solutions, while for w = 14 
it is of the order 10~^, which is extremely small and not visible on the scale of the plot. In 
Fig. ^3 we give the same information but in a diagram of current vs w. The lines correspond 
to the maximum current for each lobe and below each line there are two solutions. Thus 
the solutions ao and ar have the same maximum current (starting from zero current) and 
correspond to the right lobe of Fig. ^d, while u and al to the left one. In order to make sure 
that we get all solutions we scan over m (with a uniform and fine grid), and the current is 
obtained from (pTSl). As expected the maximum current for large w is accurately estimated 
by (|18|) down to w = 4, while for small w it varies linearly. 



B. Fixed m solutions 

Another possibility exists if w > vr, so that we can fit in the length exactly an odd 
number of half periods. This automatically satisfies the antisymmetric boundary conditions 
due to the current. Then, there exists a fixed (sometimes more than one depending on 
the length) m = m* for which w = 2K{m*). In fact every time w increases by 27r there 
is an extra solution arising. Thus for 7r(2n -\- 1) < w < 7r(2n + 3), we have solutions 
aX w = 2K{'m), ■ ■ ■ ,2{2n + l)K{'m) with ra-different values of m*. By shifting xq we can 
obtain a range of possible currents / while always satisfying the boundary conditions at zero 
magnetic field. In Fig. |]a we plot the current as a function of Xq, for w = 10, where we 
expect two solutions of this type. The corresponding values of m* are m* = 0.999272 (from 
w = 2K{m*), see curves 1/ and Ir in Fig. |^a) and m* = 0.213839 {w = GKlm*), see curves 
0/, Or in Fig. |^a), with the maximum currents being Jq ~ 4 (at xq = and Ji = 1.8 (at 
Xq = ^) correspondingly. These also coincide with the maximum currents obtained by fixing 
Xq = —K{m), K{m) and varying m. It should be remarked, though that they correspond 
to different solutions as we increase the current (/ < Imax)- This can be seen by the different 
fluxon content of these solutions in Fig. At the maximum current value though they 
coincide. In comparison, the solutions in the previous subsection correspond to zero fluxon 
content Nj = 0. 
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C. Solutions for w = 10 

We will present in more detail the calculations for w = 10 since this is a common length 
in experimental design of long junctions and one can clearly see the multiplicity of solutions. 
In Fig. P we present all the solutions (for constant Xq and constant m) for = at three 
different currents in order to follow their evolution. The branch with a half period solution 
for w = 2K{m) (see Ir and 1/ in Fig. ^-i) has a maximum current J = 4 which is the same 
as the maximum current in Fig. In fact at the maximum current (Fig. ^) we have the 
coalescence of four different solutions. In the third column of Fig. ^ (i.e. g, h, i) we see the 
four solutions being different at / = (see g) but converging to the same solution (modulo 
27r) as / approaches the 7^4 (see i) value which is the maximum current for all four 
solutions. The four solutions come in pairs: two from the pair with w = 2K{m) discussed 
earlier (i.e. Ir, 1/) and two from the right lobe of Fig. Qd, i.e. ar and aO in Fig. ^-i, discussed 
in the previous subsection. The other pair of solutions with 3 half periods when w = 6K{m), 
i.e 0/, Or in Fig. |^a-c have a maximum current near / = 1.8. For higher currents (above the 
value at c) it jumps branch and converges to the solutions of the left lobe of Fig. since the 
two pairs of solutions are quite close as can be seen from the plots in Fig. ^ and Notice 
that the currents are different for the two plots. We should also point out (to be discussed 
in the next section) that by slightly increasing the magnetic field the w = 6K{m) solutions 
show an interesting bifurcating behavior with a jump in the maximum tunneling current. 
All solutions as seen in Fig. ^ come in pairs with opposite fluxon content. Thus on the line 
Nf = (zero flux) we also have two pairs of solutions up to / = 2.4 (point A in Fig. ^b, 
where only the point at the maximum current is shown). One pair is the curves u, al in 
Fig. e, f. A second pair goes up to / = 4.0 (point B), i.e. the curves aO, ar in Fig. ||g-i. 
The solutions of each pair are different as can be seen in Fig. |^a, d, g or Fig. |]b, e, h and only 
coincide at the maximum current. Notice that u and aO are simply displaced by tt, while al 
and ar have opposite signs. With increasing current, though, they evolve very differently. 

IV. NO CURRENT CASE (/ = 0) 

In the absence of external current Eqs. (|^J^) reduce to the following 

2\fmcn ( xo + -^1"^ ) = H-, (19) 



l^^cn \ xq — ^l^^j = ^1 (20) 

which determine the parameters m and Xq that characterize the periodicity and phase 
shift of the static solutions. This is seen from Eqs. ([T9| , |20|) and the periodicities of the elliptic 
functions. Notice that these are the only three possibilities leading to three branches, if we 
consider solutions, where only m varies and Xq is fixed to Xq = or xq = 2K{m). Here we 
need Xq = 2K{m), and not Xq = —K{m) as in the zero field solution, due to the symmetric 
boundary conditions for the magnetic field. At the same time we also have solutions where 
m is fixed and xq is varied continuously with the current. 
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A. Fixed xq solutions 



We start by giving the three branches. 



Branch I: An obvious choice in (p!9| . |20| ) is Xq = 0, so that for m < 1 

1JJ 

H = cn{ — \m) , < m < 1, (21) 

cos$o(a;) = 2 m sn^(x|m) — 1, < m < 1, (22) 



sin $0(3;) = ~2 i/m sn{x\m)dn{x\m), < m < 1. (23) 

Another possibihty is xq = 2K{m) for the case that K{m) < ^ < 3K{m), or more 
generally when cn{^\m) < 0, since we are limiting ourselves to H > 0. This means that 
every time w increases by 2tt we introduce two extra solutions. To put it in other words the 
function H{m) in (^Tj) is highly oscillating for large w. 

Branch II: For Xq = and m > 1 we use the notation m = 1/m and the transformation 
rules of elliptic functions to obtain 

H = ^ dnl-^\rh] , < m < 1, (24) 
y/m \ 2 vm / 



cos<l>o(a;) = 2 sn^ ( a;^=|m ) - 1, < m < 1, (25) 



sin $0(3;) = ~2 sn ( x—^\m \ cn (x |m] , < m < 1. (26) 



Branch III: Taking into account the period and symmetry of the nd elliptic function 
we can also put in Eqs. (|l|,|OD Xq = ^^(^) with m > 1 (it cannot be satisfied for m < 1) 



to obtain 



H = 2J^^nd\ 7^|m , m < 1, (27) 



cos$o(a;) = 2 cd^ (x^\rh] - 1, (28) 



1 



sin $0(2;) = 2 y/l — m cd x—^\m sd x—j=\m . (29) 

V \/m / V \/m 



The branches II and III were obtained from Eqs. ([T9|j20|) by assuming that the modulus 
m > 1 and putting rh = ^. The expressions (^l]) and (|^) can both be written in the form 
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H = cn{ — \m) , < m < oo. (30) 

In this case branch II is described by Eq. (^) at m > 1 which reduces to (|25|) by using the 
transformation formulas (see Ref. ||^ ) when the modulus is greater than unity. Again for 
large H{m) is strongly oscillating, even though it remains always positive. 

In Fig. ^ we plot the fluxon content A'^^, at / = for three lengths, that show different 
patterns and evolution and the introduction of multiple solutions with length. In the first 
row we plot the magnetic field {H vs m) and we see that for the same H we have different 
flux. This means that they correspond to different solutions with different screening currents. 
Comparing the second and the third rows we see that the modulus m (at zero current) is 
a much better parameter than the magnetic field H to characterize the solutions, since the 
flux in this case is unique except for a symmetry multiplicity. Also notice that the curves 
should be symmetric about the horizontal at the zero level (two xo values), but the plot is 
not completed to keep the vertical scale shorter and avoid optical complexity to the eye. 
Only for branch I we show both curves (for xq = and xq = 2K). On the other hand the 
magnetic field plots exhibit some interesting changes of slope which for very long junctions 
alternatively correspond to stable and unstable regions of solutions. The change of slope is 
especially apparent for w = ^ but the alternation of stable and unstable regions also exists 
for small lengths as will be discussed later. 

Notice that at the points of slope change the fiuxon content is an integer for both branches 
II (light dashed) and III (dark dashed). Actually for stronger H the curves in Fig. ^i will 
look just like in |^c. Notice also that the symmetry in figures (b),(e),(h) will correspond to 
an antisymmetric form in the plot of flux (Nf) with H. The oscillatory form of flux vs H 
is understood by looking at the relation of H and m as plotted in (a),(d) and (g) for the 
three lengths. Notice the evolution with the creation of lobes for small m, whose number 
will increase with w as discussed. In the w = ^ case we also have extra solutions with m 
fixed which are not shown in the plot, but will be discussed in what follows. Also for higher 
w the Nf plot becomes more complex and as a particular case we discuss the w = 10 length 
in the next subsection. 

B. Magnetic flux for u; = 10 

In Fig. ^ we plot the fluxon content Nf, for a junction length w = 10, at zero current 
as a function of, the magnetic fleld H in (a) and equivalently the modulus m in (b). We 
see that the plot in (b) is essentially single valued, while the two curves correspond to the 
choices xq = and xq = 2K{m) (of branch I), which give solutions with opposite flux. The 
corresponding plot with H is quite deformed (due to the periodic relation between H and 
the modulus m). In the plots the lines are the results of the analytic solutions and the 
symbols the numerical simulation results. In the second case we have to try different runs 
to complete the curve. As we see the plot of Nf vs m can be continuously traced by varying 
the modulus m in the analytic solutions. Then one can obtain the magnetic fleld H from 
m using analytic expressions (shown in Fig. ^) and also trace the curve in (a). When, 
however you do numerical simulations, using the magnetic fleld as a varying parameter, you 
can trace only the part of the curve with the same slope. When you reach an extremum in 
the magnetic fleld (see points g, I, i, k, r, q, - ■ ■ in Fig. ||a), where the slope becomes inflnite, 
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the iteration procedure for the branch continuation with increase of H does not converge. 
Then you must decrease the value for H to trace the negative slope curve. Thus one needs 
five tracings back and forth in H to obtain all the solutions of branch /, i.e. the curves 0/, 
Or, u, 11, Ir in Fig. ||a. 

In Fig. ^ we show all three branches : branch I as given by the curves o — i — g — e for 
Xo = (solid line) and o — k — I — m) for xq = 2K{m) (dotted line); branch II is given by 
the line e — c — a — ■ ■ ■ (long dashed line) ; and branch III is given by the line o — r — q — p ■ ■ ■ 
(dashed line). Near H=0 for 1=0 we have several solutions four of which have the same 
fiuxon content Nf = 0, i.e. al, ar, aO, u, and four with different Nj at points e,m,x and 
y, i.e. 11, Ir, 01, Or. Notice that the point e is slightly to the right of the Nf axis. This 
is because, for the particular value of w = 10 near m = 1, the cn{^, 1) elliptic function 
behaves like sech(|^) so that for large w it is small and positive. 

When we increase the current / and fix if = 0, the points x, y of Fig. ||, will give us the 
curves Or and 01 in Fig. ^ for the flux. The corresponding points in the neighborhood of 
e and m will give us the curves 1/ and Ir in Fig. ^ At the point o there are two solutions 
with constant m, i.e. al, ar in Fig. ^ (which for H = are part of the curve Nf = 0). 
With increasing current one of them (i.e. al) reaches the maximum current at point A in 
Fig. with I = 2.5. The other (i.e. ar) goes to B in Fig. ^ with / !^ 4.0. These solutions 
have m = 0.88299 where cn{w\m) for the given w vanishes. They have opposite flux because 
they correspond to the two possible values of xq = —K{m),K{m). They are equivalent to 
the solutions in the middle point with / = in Fig. ^d. In the point o there are two more 
solutions from which one belongs to the stable branch o — r (i.e. aO) and the other to the 
unstable branch i — o — k (i.e. u) in Fig. |^. 

C. Fixed m solutions 

The first branch (I) has also solutions with fixed m = m* = 0.88299 if w > 27r. The 
value of m* is obtained from the condition ^ = 2K{m*), so that we have an integer number 
of periods {4K{m)) for the elliptic function in the junction length w. This automatically 
satisfies the symmetry requirement for the boundary conditions. The magnetic field is 
determined from the position phase parameter Xo (with H a periodic function of Xo with 
period AK{m*) = w). It is given by 

H = 2^/m*cn{xo + 2K{m*)\m*) 

and is presented in Fig. ^ for w = 10. The maximum value of H for this branch is H = 1.8 
at Xo = ±2K{m*). The two signs correspond to the al and ar curves. Notice that these 
solutions at zero current have zero fiuxon content {Nf = 0) over the whole extent of the 
magnetic field for which they exist. These solutions also exist for w = ^n, but are not shown 
in Fig. 01, i. For larger w we expect more pairs of solutions, i.e. a pair for each increase of 
w by 27r. Thus for An < w < Git there are two pairs of solutions. 

D. Maximum tunneling current 

In Fig. |l^ we plot the maximum tunneling current as a function of the magnetic field 
for the three different lengths. We see that for w = ^ there are two curves (like ao and 
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Mo in (a)) for the maximum current, one of which is stable and the other unstable. There 
are, however, abrupt variations of the maximum current. Thus at the end of the stable 
(0-1) branch (line a — b) there is an unstable (1-2) branch (line c-d), which however has a 
discontinuity at about H = 3.3 in I^ax (see vertical arrow). The same happens at the left 
of the stable (1,2) branch d — n — c — m where its continuation has a discontinuity of 0.4 in 
Imax, at about H = 2.6. We see that at the Imax there are two curves superimposed with the 
same Imax but they start from different solutions at / = 0. Thus one has to be very careful 
when tracing numerically the Imax vs H curve. For intermediate current values (/ < Imax 
at a fixed H) you must check to the right of c which solution branch you follow. Thus if 
we trace for the Imax the stable branch at / = from if = (i.e. ^o{x) = 0) to the right 
we follow the curve a — p — b — c — n — d while if we start from the unstable branch at 
H = (i.e. $0(3;) = vr) we follow the curve a— p — m — c — n — d---. 



The plot ioiw = f (Fig. [Tgb) is a bit more complicated. This to some extend is caused 
by the loops in Fig. |^. Thus when we trace from H = 3.0 to the left, we follow the curve 
a ~ b — c with a jump to d then d — f — e — f — g — f — h, followed by a jump to a point 
symmetric to c and from then on following a symmetric path which is not shown in the 
figure. The corresponding path for w = ^ (Fig- P^c) is a — 6 — c with a jump to d then 
e — f — g, to h and from there on a symmetric curve. Notice that in this case, the extra 
branch e — f ~ g has a lower peak current from the previous cases. 

In the above Im,ax diagrams the existence of jumps as we scan the magnetic field and 
increasing the / value (at fixed H) implies a dependence of the final solution at Imax on the 
initial condition and the path of approaching it. This becomes clear if we connect it with 
the morphology of the I and H contours in Fig. ^ for the case w = 37r/2. 

As we see there are four paths to reach the point / (notice corresponding points in both 
Figs. and pi]|b) if we fix if = while increasing the current. These paths are along the 
curves xq = 3i^(m) (which is equivalent to —K{'m)) and the vertical line m = m* = 0.84. 
From these only the one from the left of / along xq = 3K{m) (up to m = 1) is stable. This 
corresponds to the solution oq. The other three paths will give solutions u, 11, and Ir. The 
last two are along the vertical line and u along xq = 3K{m) from O to f. Notice that the 
whole vertical axis (i.e. m = 0) corresponds to the single point {H = 0, i = 0) in the I-H 
diagram. 

From the contours of H around the points e and g it is clear that at these points we have 
extrema of H which also fall on the curves xq = 4:K{m) and xo = 2K{m) where i = 0. Let 
us take another look in Fig. |TD|b. If we start from the point a (with Nf = 2) by decreasing 
if at i = we reach the point c (with A'^^ = 1) along the line Xq = (i.e. branch II) in 
Fig. 1^. The continuation of branch II through m < 1 is branch I which goes up to point e 
in Fig. ^0|d. In the rest of the curve, i.e. when m goes from e to O, the magnetic field is 
reversed from -0.6 (at e) to zero. The range in H from c to m corresponds to two different 
paths in the Xq — m diagram. From the topology it is clear that they end up in different 
maximum current as H is kept constant. The path mc (above m) has its maximum current 
to the left of the line, while the path cm (i.e. below m) has its maximum current to its 
right. Notice that in Fig. [IO|b, in going from a to c (along branch II) you cross the H value 
at m. The corresponding point in the Xo-m diagram is a different one (m') on the m-axis 
between a and c. Finally the point b can be reached from several paths. A similar analysis 
can be given in all cases, but we chose a single length as a point of illustration. 
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Below each case in Fig. |T0| we plot the corresponding flux vs H at the maximum current. 
These curves resemble the ones at zero current in Fig. ^(c, f, i). One remark is that while 
at / = the Nf vs if is a continuous curve, at Imax the flux shows discontinuities the same 
way that the maximum current was discontinuous. Comparing Figs |^f and 0e we see that 
the branches om (o being the origin) and ca are only slightly modifled. The branch though 
e — c (in Fig. |^) is folded onto e — o — d (in Fig p!0|e). The same happens for w = ^ between 
Figs. 1^ and lOf. In the plot for the flux we have not plotted all the branches. For a more 



complete plot of the branch II and III solutions see the case w = 10 in Fig |Tl] since they 
have the same approximate structure. 



In Fig. |Tl]a we show the maximum tunneling current and the flux (at Imax) as a function 
of the magnetic fleld, for w = 10. At this length we are already at the limit of long length 
branches. In this case when going from right to left we trace the Imax through the points 
a — b — c — d — e, jump to f — g — h — i — j — k — l — m with a symmetric continuation. The 
flux (Fig. pA[b) also shows a similar folding as for the case w = The letters correspond 
to the ones in Imax vs H plot. 



V. ANALYTICAL STABILITY ANALYSIS 

Next we obtain some analytic estimates for the stability regions for the zero current 
solution. The observations obtained by these estimates will verify and extend our results by 
solving numerically the stability eigenvalue problem in ([l^). Substituting Eq. (^) into Eq. 
(p!T|) we obtain the Lame eqn. 

-X" +\2m sn'^{x + xo\m)-l\x = \X. (31) 



Numerical solution of (^T]) with the boundary conditions in ([T^) gives us all the eigenvalues 
and we can check the stability of the static solution. One can gain some insight into a nec- 
essary bound for stability from the following analytic considerations, which we also compare 
with the numerical results. 

Let us remark that for three values of A = m — l,0,mwe can give an explicit analytic 



form for the corresponding solutions of (0), which are 



Xq = dn(x + Xo\m) , Xq = m — 1, (32) 
Xi = cn{x + xo\m), Ai = 0, (33) 

X2 = sn{x + xo\m) , X2 = m, (34) 

while other eigenfunctions of Eq. (|3T| ) have much more complicated forms. It is worth 
stressing that the functions in (|32| , ^|34|) cannot be called eigenfunctions of our problem 
in ([Tl])-(|T2|), because they do not satisfy the boundary conditions (p^). Taking now into 



account Eq. (|T2D we can flnd the curves of neutral stability, i.e. the critical relationship 
between parameters {w and H or I in our case) when the problem (|TI|)-(|T2|) has an eigenvalue 
A = 0. 
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In the following we shall examine the implications of the above three analytic eigenfunc- 
tions on the stability of the static solution. We must bear in mind though that it is not 
sufficient that one of the three above eigenvalues is zero or negative, but on top we must 
satisfy the boundary conditions. Even in that case we can prove instability but not the 
reverse, which can only be done by numerical solution of the eigenvalue problem. We will 
see, however that some of the conclusions will be very useful. We will examine the situation 
for each of the three branches separately. 

1st branch: In this case it is not easy to get analytical estimates. We do not have an 
extra neutral stability since m < 1 and in any case it can be considered as a continuation 
of branch II. 

Ilnd branch: In this case the three eigenfunctions of interest are 



cn 



m x\ — 
m. 



An 



m 



1 > 0, 



(35) 



Xi = dn 





= dn 


X 


\pm X — 


—j=\m 


m. 







Ai = 0, 



(36) 



Xo = sn 



mx 



1 

m. 



Ao = m > 0. 



(37) 



Again, since m > 1 only Xi is of interest. It is an eigenfunction of the linearized problem 



if 



/ w 
sn I — ^\m 



m 



(38) 



or equivalently if 



w 



2jK{m), J = 1,2, 



(39) 



Substituting (|39| ) into ( P4[) we get two families of curves of neutral stabihty, where again 
we can distinguish two cases for even (j = 2n) and odd (j = 2n + 1) values of j. For 
even j we have Nf = 2n + 1 and for odd Nf = 2n , i.e. odd and even number of fiuxons 
correspondingly. The respective values of w are given by 



w 



H 



-nK 



4 \ 



with ri = 0, 1, 



(40) 



w 



2{2n + l] 



K 



Nf = 2n + l , with n = 0, 1, • • • . 



(41) 



Ilird branch: In this case m > 1 and we can write the eigenfunctions in the following 
form ^ 

Xn = cn \/m{x + Xo)\ — 



m 



An = — 1 > 0, 
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Xi = dn 



m{x + Xo)|- 



m 



dn 



m x\ — 



m 



dn 



\m 



Ai = 0, 



(42) 



sn 



m{x + Xo)|- 



m 



Xo = m > 



In ( ^21) we used a standard transformation of elliptic functions so that their modulus is 
less than unity, and then substituted xn = \ —Ki—). We also used the notation ffi = —. 
Since Xq and A2 are always positive we only need to consider Xi. For Xi in to be 
an eigenfunction of (|lT])-(|T2|) we must have again the condition ( p8D or the equivalent (p9[). 
Here we can distinguish two cases for even (j = 2n) and odd (j = 2n + 1) values of j. It 
can easily be verified again that for even j we have Nf = 2n and for odd Nf = 2n + 1 , i.e. 
even and odd number of fluxons correspondingly. 

Substituting Eq. (^Uj) into Eq. ( P7D we obtain two families of curves of neutral stability 



in H, with the values of w given for n 



1,2, 



by 



An 



w 



1 + 



-.K 



1 + 



N = 2n 



(43) 



and 



4 / 4 \ 

w = —{2n - \)K ( — I , N = 2n. 



(44) 



Eq. (1^) is valid for any H > but (^41) only for H > 2. As we will see in the section of 
the numerical evaluation of stability the above families of curves will in fact compare very 
well giving the boundaries of stability. 



A. Numerical stability results 

To check the stability of the solutions discussed and examine the validity of the analytical 
stability results we have calculated the eigenvalue spectrum for small oscillations around the 
solution ^o{x). In Fig. |T^ we plot the four lowest eigenvalues of Eq. ( PT]) for w = 10 
and / = as a function of the parameter m for the I and II (Fig. |12|a) and the III 



branch (Fig. |T2|b) correspondingly. Stability requires all eigenvalues to be positive while an 
increasing number of negative eigenvalues denotes a higher degree of instability. 

The regions in m with all positive eigenvalues correspond to the stable fluxon solutions 
and are separated by regions in m with unstable solutions. Often stable solutions correspond 
(in the plot of Nf with H) to the branches with positive slope, line ec, etc in branch II 
of Fig. ^ or qp, etc in branch ///. This is not the case though of small w (see Fig. |^ or 
strong magnetic fields. All the solutions in branch / (m < 1) are unstable, with several 
eigenvalues being negative. In all expected the lowest mode has no nodes and is 

symmetric. It can have however several lobes, reflecting the number of fluxons that show in 
the unperturbed solution ^o{x). Also when a higher mode eigenvalue vanishes the lowest 
mode reflects this and reforms by creating more lobes, and this effect can be strong when 



eigenvalues cross each other. In comparing Figs. 12 a and b we see that the regions of stable 
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solutions in m for branch // are regions of unstable solutions for branch ///, while the 
bounding values of m are the same in both cases. This is consistent with the analytic 
expressions (26) and (28). Of course they correspond to different solutions due to different 
Xq. In branches // and /// there are at most two negative eigenvalues, while in branch / 
[m < 1 in (|34D] there is a region with four negative eigenvalues. The analytic formulas for 
stability are entirely consistent for the points in m where instability sets in. In fact in Fig. ^ 
we show the four lowest eigenvectors for m = 1.01012 where Ai = 0. The corresponding 
eigenvector is fitted with equation (25) and as expected for m ~ 1 the dn function behaves 
like a sech(y^a;). The same is true for the branch ///, where the lowest mode can be fitted 
well with (31). 

In Fig. [ij we plot in the m vs w diagram the lines from conditions (^), so that each 
line corresponds to solutions with an integer fiuxon number. Thus the range of m values 
between two lines for a given w corresponds to the {j,j + 1) branch of both II and III cases. 
In Fig. |l^ we plot the same information but in an H vs w plot for II (Fig. 15a) and III 
(Fig. p!5|b). Due to the oscillatory relation of H and m the extrema of the magnetic field at 
1 = for each branch do not exactly coincide with the stability boundary lines. This means 
that the branch ends on these lines but in the intermediate it might reach H values slightly 
outside this range. If viewed as a function of w for a fixed m (or m) then the corresponding 
magnetic field varies periodically, with a period in w equal to A\f^K{jn), i.e. for each m 



it covers the if— values between every other curve (Aj = 2). We also notice from Fig. |15 
that the width in H of the (0,1) branch is almost independent of w for large w, while for 
w ^ oo many branches tend to coalesce in the same interval of H for both II and III. One 
can understand this case by using the pendulum analogy and realize that for large w the 
solutions correspond to trajectories that pass near the separatrix points. This can also be 
seen from (0). In order to satisfy the boundary condition for w oo, the important values 
of rh are quite close to m ^ 1, since in that case K{m) oo. On the other hand for high 
magnetic fields H ^ oo the corresponding values of rh are near zero. 



VI. EXPERIMENTAL RELEVANCE 

In this section we relate some of our theoretical and numerical results to the experimental 
techniques and data. Since the behavior of the maximum tunnel current is of importance for 



junction characterization, we start by plotting in Fig. |T6| the numerically estimated (using 
the iteration procedure described in section 4.2) values of Imax for three different lengths 
(8.24, 9 and 10) and the associated experimental data extracted from figure 5.10 in [Q. The 
best fitting seems to be for L = 9.0, i.e. slightly different from L = 8.24, as was determined 
by analysis of the experimental data. The discrepancy is due to the fact that the critical 
current density is not homogeneous in the experimental sample and the analysis used in 
the experiment it would be valid for a larger length junction. An exact knowledge of the 
inhomogeneity can give a more accurate profile but this is not the purpose of this paper. 

Besides the relatively good agreement of the maximum current Imax value for each H, that 
verifies our numerical and theoretical analysis we can also make the following observations: 
The experimental data for Imax seems to try to follow the stable branches. After crossing of 
Imax lines it seems to approach some "bifurcation" points where it becomes easy to fall in 
another branch, but a careful experimental quasistatic scanning of the magnetic field and 
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current ( as was done in the numerical simulations used to obtain the displayed data), will 
enable to successfully trace the whole stable branch experimentally. Then one can pass over 
these 'fuzzy' bifurcation points and be able to select the appropriate continuation branch. 
To work within a given branch with low Imax is of interest for low energy (or current) devices. 
The quasistatic scanning will also help to elucidate the physical nature and the practical 
consequences of such bifurcation points. 

There is a close relation between the experimental and the computer simulation method- 
ologies for determining the whole Imax line of a device. One of these methodologies is 
described (as a numerical scheme) in It is based on the failure of the convergence 

of the associated iteration method when the bias current exceeds the maximum valued. We 
should mark that in this case one has to solve a large set of PDE problems associated with 
continuation points on the Imax line through fine tuning of / and H along the boundary 
line. This requires significant amounts of computer power since its is prone to the effect of 
hysterisis. A similar approach is used in experiments. In this case tracing Imax corresponds 
in configuring the device on the border line of the branch. Specifically slightly higher bias 
currents switch the device from the pair tunneling to quasiparticle mode. A similar approach 
(used here) is to configure the device so it corresponds to a point on the H-axis (1=0) and 
slightly increase the current until the above mode switch will be detected. It is possible but 
not easy to realize such "initial configuration". For both cases one has to introduce some 
initial flux conflguration. 

The main difficulty that arises from the hysteretic behavior can be avoided if one looks at 
Fig. |l^ where we replot Fig. |^ but with a rescaling of the vertical axis by K{m), so that the 
curves xq = K{m) now become horizontal lines. As we mentioned earlier the parameters xo 
and m define uniquely the fiuxon distribution of the solutions. So the shaded area enclosed 
by the curve through the points a, b, and c (with upper half / > and lower half / < 0) 
is the stable region corresponding to the (1 — 2) branch in the Imax{II) diagram (see the 
corresponding a, b and c points). It is actually separated from the other stable regions. 
One can go however from one stable region to the next, by moving quasistatically along the 
Xq = K line. 

Then it is clear that for an experiment one might select the starting configuration of his 
choice and follow an appropriate path quasistatically to another stable region and to the 
maximum current so that the whole procedure is convenient. Note that keeping H (or /) 
constant corresponds on walking on a particular contour line. 

For the procedure described above where the H is increased or decreased monotonically 
the scanning in the I—H plane suffers from strong hysteretic phenomena that are apparent in 
both the computer simulations and the actual experiments. Based on the analysis presented 
in this paper an alternative way free of such phenomena can be very naturally proposed. 
Specifically, as seen in Section 4.2 (see in particular Figs. 8a-c) one might consider searching 
for the Imax on the / — Nf plane where its is mostly singled valued. The search methodology 
will remain the same as before and therefore it can be easily done on computer simulations 
by making simple modification on the existing software. Nevertheless it is not clear how 
this can be done experimentally since it requires a manipulation of both the current and 
the external field. Another way to keep the fiux constant and this can be done by applying 
a non-constant magnetic field that varies trying to keep the distance between the fiuxons 
constant. It remains to be seen how easily this can be done in practice. 
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VII. DISCUSSION 



In the preceding sections we have presented a theoretical, numerical and experimental 
study of the various static solutions of 1-D Josephson junctions. Our basic approach was 
the use of elliptic functions to analytically express the solutions of the associated sine- 
Gordon equation. The two parameters involved in the elliptic functions (m and Xq) were 
properly selected based on the particular form of the boundary conditions. This let us obtain 
useful analytic expressions for these solutions in particular for the cases of zero magnetic 
field if = or zero current 1 = 0. Their importance lays in the fact that one can easily 
study their stability by using simple linear perturbations. This simplicity in the stability 
analysis let us exploit the role of the geometric {w) and physical parameters [H, I, Nf) 
involved. A significant outcome of our study is the fact that the module (m) of the elliptic 
functions is a good characterization parameter that greatly simplifies the general qualitative 
and quantitative pictures of the various solutions. The use of m as a characterization 
parameter also leads to more stable, accurate and efficient numerical algorithms used to 
study various aspects of Josephson junctions. 

The analysis presented above is particularly useful to understand the sometimes compli- 
cated behavior when we try to follow numerically the different branches. In fact, this was 
one of the motivations behind this work. We are currently building a software engine to 
numerically simulate 2 dimensional window Josephson junctions of various types and con- 
figurations. Preliminary numerical experiments clearly show that although our simulation 
engine is build on top of powerful numerical continuation methods using state-of-the- 



art PDE software [|Tl|;0 very often fails if we do not fully utilize the results obtained in the 



current study. Some the most common, annoying problems encountered (even in the case 
where no defects are present) are the following: Considering bifurcation points as regular 
point and vice versa, missing bifurcation points, viewing certain turning or bifurcation points 
as limiting points and improper branch switching (e.g. 27r jumps). The present study gives 
several hints to help us drive our simulation engine with no such problems. 

Our original goal was the study of the influence of the critical current density ( Jc Jd^)) 
inhomogeneities on the tunneling current Imax- Two observations, however, made necessary 
to study the perfect junction: 

(a) We noticed that when studying window junctions, even for zero magnetic fleld, the 

maximum current starting with different initial conditions, was not always the same 
i.e. at Imax = 4.0. Several times it stopped at lower values. 

(b) Both numerical and experimental results show strong hysteresis phenomena with jumps 

between different branches when varying the external magnetic fleld. Related to this, 
the question arises whether there exists a way of analytic continuation between different 
branches? Or in more physical terms whether there is a physical parameter (in the 
place of H) whose smooth variation shows no hysteresis in Imax- 

With respect to point (a), it is clear now that the second value belongs to one of the 
unstable branches we discussed for w = 10. If, however, there are defects in the junction the 
unstable solutions might also become stable and therefore are of interest [Q. Another way to 
stabilize solutions is by high frequency fluctuations (of small amplitude) in a way similar to 
the Kapitza inverse pendulum problem ||15|. This can also be achieved by small wavelength 
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spatial variation of the critical current density . For remark (b) in the undefected 1 — d 
junction one has the advantage that the analytic solution is known and the choice of m and 
Xo, pin uniquely the proper solution. Thus one can follow, smoothly the solution if we look 
at Imax as a function of the magnetic flux. This way we can avoid hysterisis by choosing a 
proper initial condition at / = and increase / to Imax- 

If, however, one uses the magnetic field as an input parameter, strong hysteretic phe- 
nomena are observed, due to the non uniqueness of the relation between H and m for large 
w. The multiple solutions (for fixed Xo due to symmetry) correspond to different fiuxon 
content. This non-uniqueness will disappear for large H where in fact the junction behaves 
as if it is a short one and you recover the diffraction like pattern. Also an increase of the 
temperature makes the junction to behave as a short one, with non-overlapping branches. 
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TABLES 



TABLE L Different types of curves due to various selections of m and xq. 
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FIG. 1. Plot of the maximum tunneling current as a function of applied magnetic field for a 
short junction with: (a) w = 1.0 and (b) w = 10. 
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FIG. 2. Constant H (a) and / (b) contours in the (m, xq) plane. The curves with the symbol 
O are for = in (a) and / = in (b). The signs " + " or " — " give the sign of / and H. The 
letter symbols signify the same points in the / and R contours. 
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FIG. 3. (a) Plot of I{m) (using Eq. (13)) with xq = —K{m) or K{m) and three values of w. 
(i) w = 27r/3, (ii) w = 37r/2 and (iii) w = 57r/2. (b) Same as (a) with u; = 10 (dashed line). The 
right lobe in (b) is also shown in an expanded scale for m (solid line). The small m scale is shown 
at the top of the figure. 
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FIG. 4. (a) A diagram of the different solutions in the space of parameter w and m with 
xq = zizK{m) for the solutions in equations (13)-(15). (b) The same as (a) but presented in the 
space of w and I (instead of m) . 
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FIG. 5. (a) Plot of I{xo) ior H = 0, w = 10 and fixed m with (i) m* = 0.999272 (curves 
11 and Ir) from w = 2K{m*), (ii) m* = 0.213839 (curves 0/ and Or) from w = 6K{m*). (b) 
Fluxon content Nf as a function of the bias current I for the different solutions with fixed m: (i) 
m=0. 999272 with one half period of the elliptic function (Ir and 11). (ii) m=0. 213839 with three 
half periods of the elliptic function (0/ and Or). 
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FIG. 6. Plot of the phase distribution for all the solutions at if = and three values of 
the current (different for each line), (a, d, g) I = 0, (6, e, h) I = Imax/'^-, (c, /, i) I = Imax 
(where Imax is different for each case). The three columns correspond to: (a, b, c) the two solutions 
(0/, Or) with m* = 0.213839. {d, e, f) the two solutions (n, al) from the left lobe of Fig. 3b and 
{g, h, i) the two solutions (Ir, 1/) for m* = 0.999272 and the two solutions from the left lobe of 
Fig 5a (i.e. aO, ar) 



27 



w=27i;/3 



w=37c/2 



w=57t/2 



4.0 - a 



---^^1 4.0 - d 



S 2.0 

n: 

0.0 

-2.0 
0.0 



1.0 - 



l,x =2K 



- 20 - /r 

- 0.0 



2°m 



-2.0 



-1.0 



0.0 



2.0 - e 



0.0 



0.0 c 



-^-2.0 
4.0 0.0 



2-°m 4° 




- 3.0 



- 1.0 



1/ 



-^-1.0 
4.0 0.0 



-1.0 



-2.0 0.0 2.0 4.0 
H 



3.0 - f 
-1.0 




-1.0 1.0 3.0 5.0 -2.0 0.0 2.0 4.0 
H H 



FIG. 7. Plots oi H Ys m (top figures) and of the fluxon content Nf for I = as a function of 

3 ' 2 ' 2 



m (middle figures) and H (bottom figures), for a junction of length w = 
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FIG. 8. Plot of the fluxon content Nf for / = and w = 10 as a function of (a) H and (b) m. 
In (c) we give H{m). 
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FIG. 9. Magnetic field H{xo) at / = for the constant m = 0.882992 solution. 
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FIG. 10. Maximum tunneling current as a function of H for three different lengths (a) w = 
(h) w = (c) w = ^ and A^/ vs H at the maximum current (d),(e),(f) correspondingly for the 
three lengths. 
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FIG. 12. Plot of four lowest eigenvalues as a function of m for (a) branches / and // and (b) 
branch ///. 
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FIG. 16. Numerical results for Imax vs H for L = 8.24, L = 9 and L = 10 and the experimental 
data from [1]. In the experimental data the dots are the measured points and the dashed lines 
should only be considered guide for the eye. 
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